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ABSTRACT 

A two-dimensional (2-D) computer code was developed for 
modeling enclosed volumes of gas with oscillating boundaries, 
such as Stirling machine components. An existing 2-D 
incompressible flow computer code, CAST, was used as the starting 
point for the project. CAST was modified to use the compressible 
non-acoustic Navier-Stokes equations to model an enclosed 
volume including an oscillating piston. The devices modeled have 
low Mach numbers and are sufficiently small that the time required 
for acoustics to propagate across them is negligible. Therefore, 
acoustics were excluded to enable more time efficient computation. 
Background information about the project is presented. The 
compressible non-acoustic flow' assumptions are discussed. The 
governing equations used in the model are presented in transport 
equation format. A brief description is given of the numerical 
methods used. Comparisons of code predictions with experimental 
data are then discussed. 

INTRODUCTION 

Background 

NASA Glenn Research Center (GRC) has been involved in 
development of Stirling engines for ~25 years. GRC began 
managing the Stirling Automotive Development Program for the 
Department of Energy (DOE) in — 1 977. The DOE/NASA 
contractors were first Ford Motor Co. (for one year) and then 
Mechanical Technology, Inc. (MTI) of Albany, NY. These Stirling 
automotive engines used hydrogen as the working fluid and were 
required to operate over demanding automotive driving cycles. 
This work continued into the early ‘90's (Ernst and Shaltens [1]). 
Engines were demonstrated in automobiles and trucks and future 
development looked promising. However, due to decreases in oil 
and gasoline prices and improvements in spark ignition engines, 
the automotive manufacturers did not choose to develop the 
Stirling technology for the automotive market. 

Development of Stirling engines for generation of space 
auxiliary power began in -1983 w F ith MTI as the primary 
contractor. These engines used helium as working fluid, were 
25 kWe (12.5 kWe/cylinder) designs with electrical power 
produced by linear alternators and used gas bearings. Back-to-back 
cylinders with synchronized pistons damped engine vibrations. 


One of the problems encountered in Stirling engine 
development by GRC and MTI w f as the limited accuracy of the 1-D 
engine design computer codes. Performance projections for new r 
designs w r ere usually optimistic. The last major design of a 25 kWe 
(12.5 kWe cylinder) Stirling engine for space power was the 
Component Test Powder Converter or CTPC (Dahr [2]). Because of 
previous power shortfalls wdth other new engine designs, MTI and 
NASA chose to design the CTPC with a 20% margin on pow r er. 
When tested the engine produced slightly in excess of the 25 kWe 
power goal, meaning that it produced almost 20% less power than 
predicted by the computer code that was used to design it. 

During this period of space power Stirling engine 

development at MTI (-1983 to 1993). concerns about 1-D design 
code accuracy led to fundamental research. A “loss understanding” 
program w r as started consisting of grants and small contracts for 
investigation of thermodynamic losses (Tew and Geng [3]). The 
participants were hosted by GRC at yearly Stirling Loss 
Understanding Workshops for several years. The long range goal 
w^as to improve the accuracy of 1-D Stirling engine design codes. 

The SP-100 (or Space Power- 1 00 kWe) reactor space power 
system program, for which these 25 kWe space power Stirling 
designs were intended, ended in 1993. Along wdth it. funding for 
the Stirling loss understanding program ended also. 

From -1993 to 1999, NASA support of Stirling engine 
development continued at a low level via NASA SBIR funding. 
Most of this funding was for continued development of free-piston 
engines at the Stirling Technology Company (STC) of Kennew ick, 
WA based on flexural bearing technology. GRC and Cleveland 
State University (CSU) also continued low -level development of a 
multi-D Stirling code via the thesis project reported here. 

Now r Stirling power convertors are being developed by DOE, 
NASA and STC as a substitute for the less efficient radioisotope 
thermoelectric generators for deep-space missions (Thieme, Qui, and 
White [4]). The success of STC in accumulating years of reliable 
operation on its free-piston flexural bearing machines (over 6 years 
on a 10 kWe convertor with no maintenance and no degradation in 
performance) has been a key factor. In late 1999. a 55 We STC free- 
piston Stirling was subjected to a number of tests and evaluations 
to determine its suitability for a deep-space mission (For example, it 
passed a simulated launch vibration test, while operating). As a 
result, in 2000, the Stirling radioisotope power system, or RPS 
(Furlong and Shaltens [5]), was baselined as the advanced powder 
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system that may be used on NASA missions such as Europa and 
Solar Probe missions. It now appears that long-duration Martian 
rover missions will be the first application of Stirling RPS. 

Stirling One-Dimensional Flow Design Codes 

Two one-dimensional-flow Stirling machine codes have been 
in primary use by NASA's Stirling engine contractors in recent 
years. These are the HFAST code developed by Mechanical 
Technology, Inc. (Huang [6]) and the GLIMPS code (now Sage, 
(Gedeon [7])) developed by Gedeon Associates. HFAST was used 
by MTI in design of the large (25 kWe) Stirling convertors 
discussed earlier. The GLIMPS code and now its successor, Sage, 
have been the primary 7 design codes used by the Stirling 
Technology Company. Sage is also used by Sunpower, Inc. While 
these codes have been used to develop excellent engine designs, 1- 
D codes assume uniform axial flow' and are thus believed to be 
deficient in modeling those interfaces in Stirling engines where 
significant changes in area take place. Past comparisons of code 
predictions, such as the GLIMPS/HFAST study done by Geng and 
Tew [8], show'ed rough overall performance agreement but 
differences in predictions of individual losses. It would be of 
significant value to be able to accurately characterize all major 
losses so that the design process could do an adequate job of 
trading off the various losses against each other. 

Objective of the Work 

The immediate purpose of the 2-D code development reported 
here was to develop a time-efficient compressible code for study and 
definition ofStirling machine type cylinder heat transfer/power (or 
hysteresis) losses; see Tew' [9] for more complete documentation. A 
long range purpose was to provide the basis for development ofa 2- 
D model ofa complete Stirling engine that could be used to study 
interaction of various thermodynamic losses. Even longer range, if 
computer so ft ware/ hardware time efficiencies are increased 
sufficiently, it may be possible touse a 2-D code (perhaps even a 3-D 
code) for engine design. The previous IECEC paper by Makhkamov 
and Ingham[10] suggests such design use of a 2-D code is already 
underway at the Laboratory for Stirling Engines at the Physical- 
Technical Institute in Tashkent. 

NOMENCLATURE 

Latin letters 

Cp specific heat at constant pressure 

d denotes differential of variable 

0^ hydraulic diameter (=4 x wetted area; wetted perimeter) 

D viscous force vector 

f mass force vector (due to gravity, for example) 
p static pressure 

P mean spatial pressure 

P a amplitude of mean spatial pressure 

P arithmetic mean of max. and min. mean spatial pressures 
q m entropy source due to non-zero mass sources 

Q density of continuously distributed heat sources 
R gas constant 

s entropy per unit mass 

t time 

T temperature 

stress tensor 

u velocity vector (= I U + ]V) 


U axial, x-direction, velocity 

V radial, r-direction, velocity or volume 

V Q arithmetic mean of maximum and minimum volumes 
>v = %pu rate of momentum change because of mass sources 
WJ non-dimensional work or hysteresis loss 

Greek letters 

y ratio of specific heats of fluid 
<p transport quantity per unit mass 

A molecular thermal conductivity 

p absolute viscosity 

p density, fluid mass per unit volume 

co angular velocity (rad/sec) 

£ mass source strength per unit mass 


COMPRESSIBLE NON-ACOUSTIC FLOW 


Assumptions 

Fedorchenko [11] discusses a number of subsonic initial- 
boundary-value problems that cannot be solved using the classical 
theory of incompressible fluid motion, which involves the equation 

V • li = 0 (where u - iU + jV and the unit vector i is the x- 
direction for both cartesian and axisymmetric coordinates; the unit 
vector j is the y-direction for cartesian coordinates and in the r- 
direction for axisymmetric coordinates.) Among these problems are: 
(1) flows in a closed volume initiated by blowing or suction 
through permeable walls, (2) flows in a closed volume initiated by 
a moving boundary such as in a piston-cylinder problem, (3) flow's 
with continuously distributed mass sources, and (4) viscous flows 
with substantial heat fluxes. Fedorchenko notes that application of 
the most general theory of compressible fluid flow may not be best 
in such cases because of the difficulties in accurately resolving 
complex acoustic phenomena and in assigning proper boundary 
conditions, 

Fedorchenko proposes a non-local mathematical model where 

V • u * 0 in general, for simulation of unsteady subsonic flows in a 
bounded domain with continuously distributed sources of mass, 
momentum and entropy, also taking into account the effects of 
viscosity and conductivity when necessary. The exclusion of 
sound waves is one of the most important features of the model. 

The most general form of Fedorchenko's compressible non- 
acoustic system of equations for simulation of unsteady subsonic, 
heat conducting viscous flows are the following forms of the 
momentum, continuity, energy and state equations: 

d( pu) __ _ 

■ + V • puu + Vp = pf + w + D ( 1 ) 


dt 

dp 

— + V ■ pu = & 
dt 

ds R 

1 -U Vs = — [v ■ AVT + Q ] +q r 

dt P 

f(s, P,p) = 0 


(2) 


(3) 

(4) 


When Vju = 0 (as it is for constant viscosity problems) then 


_ 1 

D - pAu + -pV(V ■ u) (5) 

3 


Note that the energy equation is written in terms of entropy, s, 
rather than in terms of internal energy, enthalpy, or temperature. 
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The simplifications in the Navier-Stokes equations used to 
eliminate acoustics and arrive at the above system of equations 
were: (1) Pressure at any time, t, and position, r, is split into a 
mean spatial pressure level that varies only with time and a 
A(pressure) relative to a reference position that varies with 
position and time: 

p(rj) = Pit) + Ap(rj) (6) 

(2) The pressure appearing in the equation of state is the mean 
spatial pressure that varies only with time. Therefore, from the ideal 
gas equation of state 

Pit) 

p- (7) 

RT(rj) 

So density is a function of mean spatial pressure level and the 
temperature field, (and is independent of the spatial pressure drop). 

Simplification of Equations for the Stirling Problem 

For the Stirling piston-cylinder problem of interest here, there 
are no mass or heat sources distributed within the cylinder volume 
and the gravity force on the gas is not of interest. Therefore the 
variables defined above: /, w. Q and q m are all zero. Also the 

ideal gas equation of state is sufficiently accurate for the helium gas 
that is used in most current Stirling engines of interest to NASA. 
Therefore the equations (l)-(4) reduce to the following set: 

dipu) „ 

+ V • pita + Vp — D (8) 

dt 

dp 

— + V • pu = 0 (9) 

dr 

ds R r 

— + U V.v = -|V AVTl (10) 

dt P 

p 

P- — (11) 

RT 

CAST AND MODIFIED CAST CODE 

The incompressible flow CAST (Computer ^\ided Simulation of 
Turbulent Flow) code, as originally received by Cleveland State 
University and NASA is documented in Peric and Schuerer [12]. 
Changes made to CAST in developing the modified CAST 
compressible non-acoustic code are documented in Tew r [9]. 

Transport Equation Format 

The solution technique used in the CAST code is based on a 
transport equation formulation of the governing equations. A 
coordinate-free form of the general transport equation is: 

^ + v(p« 0 -r^) = (12) 

Time rate 

of storage convection diffusion source 

The text under equation (12) categorizes the terms in the equation. 

Transport quantities, exchange coefficients and source terms for 
the continuity, momentum and energy equations used in Modified 
CAST are shown in Table I. Note that the energy equation is in 
enthalpy format, rather than the entropy format used by 
Fedorchenko. The turbulent kinetic energy and dissipation rate 
equations are not shown here due to space limitations. They are 
defined in Tew [9]. 


Table I: General Transport Equation — Transported Quantities, Exchange 
Coefficients, and Source Terms for Continuity, 

Energy and Momentum Equations 


Equation 

Name 

Transport 

Quantity/vol. 

Exchange 

Coefficient 

Source 

Term 

Continuity 

p (mass/vol.) 

0 

6 

Momentum 

pli ( momentum/vol ) 

P 

V ■ t st.s 

Energy 

ph ( enthalpy'.' vol . ) 

_ P_ 

clP 



c p Pr 

dt 


Numerical (Finite-Volume) Methods 

The numerical methods in modified CAST are almost the same 
as those used in CAST (Peric and Scheuerer [12]). An exception is 
the use of Leibniz's rule (Ferziger and Peric [13]) to account for the 
effect of the moving boundary (i.e., piston motion). Modifications to 
the governing equations themselves, to account for the change from 
incompressible to compressible non-acoustic flow', are discussed in 
Tew [9]. These included accounting for the non-zero V • u terms, 
adding the time derivative of the mean spatial pressure to the source 
terms of the energy equation, and using mean spatial pressure in the 
ideal gas equation so that the density field is a function of the 
temperature field and the mean spatial pressure. Spatial pressure 
drop still appears in the momentum equations to help determine the 
velocity field. Since the density is not affected by spatial pressure, 
the incompressible SIMPLE algorithm still applies. 

The governing equations are solved with a conservative finite- 
volume method (see Patankar [14]). The basic principle is: 

( 1 ) Discretize the solution domain by subdividing it into small 
axisymmetric (or rectangular) control volumes. Locate the numerical 
grid points in the center of the control volumes. 

(2) Discretize the transport equations. 

Discretization is done by formally integrating the single terms 
in the equations over a control volume. Application of Leibniz's 
rule and Gauss’s theorem yields an integro-differential equation 
relating the net increase in the transported quantity per unit time to 
the convective and diffusive fluxes across the control volume 
boundaries and to the source (or sink) terms within the control 
volume. This practice leads to a conservative method because 
boundary-fluxes leaving one control volume through its right 
boundary enter the neighboring control volume through its left 
boundary'. Since this principle applies to all control volume faces, 
the scheme becomes overall conservative. The approach is 
described in Peric and Schuerer [12] and Tew r [9] in more detail. 

AVAILABLE DATA FOR CAST COMPARISONS 

Rectenwald Computations and Kornhauser 
Experimental Data 

Recktenwald [15] computed heat transfer between the walls 
and gas inside the cylinder of a reciprocating compressor. A 
compressor cylinder contains intake and discharge valves, unlike 
the cylinders of a Stirling engine. Recktenwald used 2-D, unsteady, 
compressible equations (acoustics included) to simulate the 
compressor. In order to validate his computer code, he simulated a 
gas spring for comparison whth data generated by Kornhauser and 
Smith [16, 17]. The comparison between data and experiment w r as 
based solely on experimental and simulated values of non- 
dimensional hysteresis loss for the gas spring over an operating 
range. It should be noted that a gas spring is essentially a piston- 
cylinder which has no flou r to or from the enclosed volume of the 
cylinder. Kornhauser [18] reported further details of these gas 
spring experiments. He also reported on tests made w'ith a 
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modification of the gas spring test rig, to include a heat exchanger 
mounted on top of the cylinder, such that flow could continuously 
pass between the cylinder and the heat exchanger as the piston 
expanded and compressed the gas. This ‘Two-space test rig 
operated more like a Stirling machine cylinder than either a gas 
spring or a compressor. 

Komhauser’s gas spring test data was also used as a basis for 
validation of the 2-D Modified CAST compressible non-acoustic 
code. Hysteresis losses computed by CAST were compared with the 
experimental values of Komhauser and the calculated values of 
Recktenwald, over a range of gas spring operation. Also, since 
Recktenwald published plots of his calculated velocity vectors and 
temperature contours within the gas spring for two operating 
points, these were compared with similar plots of CAST calculated 
values for one of the operating points. Thus it was possible to 
compare the calculations of a compressible model, W'hich did 
account for acoustics in the computations against the compressible 
non-acoustic calculations of the Modified CAST code. Comparison 
of CAST and Reckenwald’s temperature contours over the cycle for 
a 10 RPM gas spring showed excellent agreement. Agreement 
between velocity vector plots also appeared to be very good, 
although due to some difference in the way the plots were made 
only a qualitative comparison could be made. These 2-D 
comparisons are shown in Tew [9]. 

Gas Spring and Two-Space Test Rig Dimensions 

The dimensions of Komhauser's test rigs are shown in Tables 
I and II. Due to a limitation of the CAST code in simulating complex 
geometries, it was necessary to make a slight change in the annular 
heat exchanger geometry simulated for the two-space test rig. The 
annular heat exchanger, physically mounted on top of the cylinder, 
w'as moved slightly outward so that the outer w 7 all coincided with 
the cylinder wall, and the heat exchanger volume was maintained 
the same in order to maintain the same volume ratio. Physical and 
simulated dimensions are shown in Table III (This small change in 
geometry likely tended to reduce the difference between test and 
data results, to be discussed below). 


Table II: Gas Spring Dimensions 


Phys i ca l Qu anti tv 

Symbol 

Value 

Cylinder Bore (Diam.) 

6 

50.80 mm (2 in.) 

Piston Stroke 

s 

76.2 mm (3 in.) 

Volume Ratio 

r v 

2.0 


Table 111: Two-Space Test Rig Dim. (Physical and Simulated) 


Physical Quantity 

Physical Value 

Simulation Value 

Cylinder Bore 

50.80 mm (2 in) 

50.80 mm (2 in) 

Piston Stroke 

76.20 mm (3 in) 

76.20 mm (3 in) 

Volume Ratio 

2.0 

2.0 

Annulus O.D. 

44.5 mm (1.75 in) 

50.80 mm (2 in) 

Annulus I.D. 

39.4 mm ( 1 .55 in) 

46.4 mm (1.83 in) 

Annulus Gap 

2.5 mm (0.10 in) 

2.2 mm (0.09 in) 

Annulus Length 

445 mm (17.5 in) 

445 mm (17.5 in) 

Min Pi st/Head Clr. 

2.9 mm (0.1 1 in) 

2.9 mm (0.1 1 in) 


Both test rigs used helium gas and the walls of the cylinder 
and heat exchanger were at a constant temperature of approximately 


294 K. 


RESULTS: COMPUTATIONS VS. DATA 


Gas Spring Hysteresis Losses 

For a gas spring, the hysteresis loss is the work that is 
dissipated by the spring per cycle at steady operating conditions; 
it's also equal to the heat generated in and transferred out of the 
spring. A good way to compare computational and measured 
hysteresis losses is via plots of dimensionless work as a function of 
oscillating flow 7 Peclet number. Dimensionless w f ork and 
oscillating flow Peclet numbers are defined, respectively, as 
follows: 


w bss=- 


jPJV 


P V 

r O v O 


O 

VY 


\ P o ) 


7-1 


V 7 


(13) 


Pe 


CO 


Po^_ D h 

4A 


(14) 


Recktenw'ald had previously plotted his calculated 
dimensionless losses on a plot of Kornhauser's data. The Modified 
CAST dimensionless losses were superimposed on this plot and 
the result is shown in Figure 1. This figure shows dimensionless 
loss as a function of oscillating flow 7 Peclet number. The CAST loss 
values were plotted at ‘"uncon'ccled" Peclet numbers as the 
obviously handwritten “X" symbols. If plotted at the corrected 
values of the Peclet numbers they would be shifted slightly to the 
right and would fall on Recktenwald's values (the solid 
diamonds). 

Five of the six pairs (of CAST and Recktenw 7 ald values) of 
calculated dimensionless losses agree well w ith Komhauser' s data. 
The one pairofcalculated points that did not agree well with the data 
is the pair shown at the highest Peclet number and corresponds to a 
1 000 RPM, 1 465 kPa gas spring operating condition. Recktenw'ald 
[15] discusses several plausible explanations for disagreement at 
high Pe^- 



0.1 1 10 too 1000 10000 

Figure 1: Modified CAST Dimensionless Losses Superimposed on 
Recktenwald [15] Plot of Recktenwald's Calculated Losses and 
Kornhauser's [18] Experimentally Derived Losses 
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Also in the mid-range of the data in the vicinity of Pe 0) =]Q y 

where there appears to be a higher and a lower data curve, the data 
agrees with the lower data curve. Komhauser [18] found that the 
higher and lower data “curves” were related to differences between 
data taken at high-pressure/low-speed (higher curve) and that 
taken at low-pressure/high speed (lower curve) which had 
oscillating flow Peclet numbers in the same range. Thus he 
concluded there was some other dimensionless parameter needed to 
resolve the data in this range. His experiment also showed that 
adding fins within the clearance volume of the gas spring 
suppressed the difference in losses of the two types of data (high- 
pressure/low-speed and low-pressure/high-speed). See Komhauser 
[18] for more discussion. 

Comparison of modified CAST and experimental P-V diagrams 
arc shown in Figures 2 and 3 for ~49 and 496 RPM respectively. 
Agreement is very good. When plotted in Microsoft Excel the 
curves were smooth. The ‘jagged" nature of the curves appeared as 
an artifact of “pasting” the Excel plots into Microsoft Word. 



Volume, cm A 3 

Figure 2: Comparison of Experimental and CAST P-V Diagrams. 
CAST used 18x12 grids, 120 time steps/cycle 
(Operating Conditions: 48.6 RPM, Mean Pressure = 555.7 kPa 
(80.6 psia), Wall Temp. = 294 K) 



Figure 3: Comparison of Experimental and CAST P-V Diagrams. 
CAST used 18x12 grids, 120 time steps/cycle 
(Operating Conditions: 495.8 RPM, Mean Pressure=223.3 kPa 
(32.4 psia), Wall Temp. = 294 K) 


Two-Space Test Rig Data and Calculation Comparisons 

Figures 4 and 5 show experimental heat exchanger heat transfer 
per unit area and annulus center-to-wall temperature difference for 
Kornhauser's [18] two-space rig. Figures 6 and 7 show' the 
corresponding Modified CAST calculated results for comparison 
with Figs. 4 and 5. The experimental and calculated heat transfer 
have different signs due to differences in definition of the positive 
heat transfer direction. 



Figure 4: Heat Transfer per Unit Area at Various Positions in Heat 

Exchanger Relative to Entrance to Cylinder. Komhauser [18] 

Exp. Data: Run #12071539, 201.7 RPM, 

1 .008 MPa Mean Pressure 

Figures 4 and 6 show' that peak experimental heat transfer unit area 
near the entrance is 4000 to 5000 W/rrf less than the 
corresponding calculated value; near the end of the heat exchanger 
the peak experimental value is about 9000 W/m" less than the 
corresponding calculated value. Thus even though the qualitative 
variations in heat transfer look very similar in the experimental and 
calculated plots, the quantitative agreement is not very' good. 

Comparison of the annulus center to wall temperature 
differences in Figs. 5 and 7 show- that the calculated temperature 
differences are smaller that the experimental values. This is 
consistent with the calculated heat transfers/area being larger than 
the experimental values. 

Tw'o-dimensional plots of the calculated temperatures, 
velocities, pressures, etc. are given in Tew' [9]. However, there are 
no experimental values available for comparison. 



Figure 5: Temperature Difference from Heat Exchanger Center 
to Wall. Komhauser [18] 

Exp. Data: Run #12071539, 201.7 RPM, 

1.008 MPa Mean Pressure 
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Crank Angle from BDC, degree 

Figure 6: Modified CAST Calculation of Heat Transfer per Unit Area 
at Various Positions In Heat Exchanger Relative to Entrance to 
Cylinder (34x20 grids, 120 time steps/cycle). For Comparison to 
Komhauser [18] Exp. Data: Run #12071539, 201.7 RPM, 

1 .008 MPa Mean Pressure 



Figure 7: Modified CAST Calculations of Temperature Difference 
from Heat Exchanger Center to Wall (34x20 grids, 120 time 
steps/cycle). For Comparison with Komhauser [18] Exp. Data: 

Run #12071539, 201.7 RPM, 1.008 MPa Mean Pressure 

CONCLUDING REMARKS 

The Modified CAST compressible non-acoustic model agreed 
well with 10 RPM gas spring hysteresis and P-V diagram data, and 
also with compressible acoustic calculations of two-dimensional 
velocities and temperatures. However, modified CAST calculations 
deviated from experimental values of two-space (piston/cylinder- 
heat exchanger) test rig data, although trends were predicted well. 

Recent preliminary comparisons of the commercial CFD- 
ACE+ code compressible acoustic calculations with CAST, at CSU 
(Ibrahim, et. al. [19]), for the two-space rig suggests temperature 
agreement within about 2% for somewhat different grids. Therefore, 
data and 2-D computations do not agree, for some reason, for the 
two-space test rig. The reasons will be explored further. New data 
will likely be required for further multi-D Stirling code validation 
efforts. 

When CAST and CFD-ACE+ calculations were made wdth the 
same type of grid and time step size, CFD-ACE+ was about 8% 
faster than CAST. Although not demonstrated here in comparisons 
against the highly developed commercial CFD-ACE+ code, 
Fedorchenko’s compressible non-acoustic technique, may have the 
capability for reductions in simulation time for those transient 
situations where compressibility must be simulated but acoustics 
are not important. 
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